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ABSTRACT 



\ Context. As accretion in protoplanetary disks is enabled by turbulent viscosity, the border between active and inactive 

(dead) zones constitutes a location where there is an abrupt change in the accretion flow. The gas accumulation that en- 
sues triggers the Rossby wave instability, which in turn saturates into anticyclonic vortices. It has been suggested that the 
Q.^ 1 trapping of solids within them leads to a burst of planet formation on very short timescales. 

^T^ Aims. We study in the formation and evolution of the vortices in greater detail, focusing on the implications for the dy- 

• . namics of embedded solid particles and planet formation. 

■ Methods. We performed two-dimensional global simulations of the dynamics of gas and solids in a non-magnetized thin 
protoplanetary disk with the Pencil code. We used multiple particle species of radius 1, 10, 30, and 100 cm. We computed 
the particles' gravitational interaction by a particle-mesh method, translating the particles' number density into surface 
density and computing the corresponding self-gravitational potential via fast Fourier transforms. The dead zone is mod- 

■ eled as a region of low viscosity. Adiabatic and locally isothermal equations of state are used. 
Results. The Rossby wave instability is triggered under a variety of conditions, thus making vortex formation a robust 
process. Inside the vortices, fast accumulation of solids occurs and the particles collapse into objects of planetary mass 
on timescales as short as five orbits. Because the drag force is size-dependent, aerodynamical sorting ensues within the 
vortical motion, and the first bound structures formed are composed primarily of similarly-sized particles. In addition 

■ to erosion due to ram pressure, we identify gas tides from the massive vortices as a disrupting agent of formed proto- 
00 \ planetary embryos. We find evidence that the backreaction of the drag force from the particles onto the gas modifies the 
£f*) , evolution of the Rossby wave instability, with vortices being launched only at later times if this term is excluded from 
\& ■ the momentum equation. Even though the gas is not initially gravitationally unstable, the vortices can grow to Q w 1 in 

locally isothermal runs, which halts the inverse cascade of energy towards smaller wavenumbers. As a result, vortices in 
models without self-gravity tend to rapidly merge towards a m=2 or m-\ mode, while models with self-gravity retain 
dominant higher order modes (m=4 or m=3) for longer times. Non-self gravitating disks thus show fewer and stronger 
vortices. We also estimate the collisional velocity history of the particles that compose the most massive embryo by the 
end of the simulation, finding that the vast majority of them never experienced a collision with another particle at speeds 
faster than lms" 1 . This result lends further support to previous studies showing that vortices provide a favorable envi- 
ronment for planet formation. 

^ . Key words. Keywords should be given 

??; 

1. Introduction poor sticking properties and a low threshold velocity for 

^ .„ , , , , ., ,. , , , . fragmentation (Chokshi et al. 1993), such that collisions be- 

The ill fate of the building blocks of planets in gaseous tween them usuall lead to destruction rather than growth 

disks around young stars stands as one of the major (Benz 2000; Siron ^ 20 04; Ormel & Cuzzi 2007). Such prob- 
unsolved problems m the theory of planet formation. lemg severd hinder th fo km . size , coagulation 
Beginning with micron-sized interstellar dust grains, co- (Brauer et al 2008a) 
agulation models predict growth to centimeter and me- 
ter size (Weidenschilling 1980; Dominik et al. 2007) in the In view of these problems, other routes for breaching 
denser environments of a circumstellar disk. Such bod- the meter size barner have been pursued. A distinct al- 
ies, however, are large enough to have already decou- ternative is gravitational instability of the layer of solids 
pled slightly from the sub-Keplerian gas, yet still small (Safronov 1969, Lyttleton 1972; Goldreich & Ward 1973; 
enough to be subject to a significant gas drag. The result- Youdin & Shu 2002). When the dust aggregates had grown 
ing headwind drains their angular momentum, leading to centimeter and meter size the gas drag is reduced and 
them into spiral trajectories towards the star, on timescales the sollds are pushed to the midplane of the disk due 
as short as a hundred years at 1AU (Weidenschilling to the stellar gravity. Although such bodies do not have 
1977a). Another acute problem is that such bodies have enough mass to attract each other individually, sedimen- 
tation increases the solids-to-gas ratio by orders of mag- 
Send offprint requests to: wlyra@astro.uu.se nitude when compared to the interstellar value of 10~ 2 . 
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It was then hypothesized (Safronov 1969) that due to the 
high densities of this midplane layer, the solids could col- 
lectively achieve critical number density and undergo di- 
rect gravitational collapse. Such a scenario has the advan- 
tage of occurring on very rapid timescales, thus avoiding 
the radial drift barrier. 

This picture is nonetheless simplistic, in the view that 
even low levels of turbulence in the disk preclude the mid- 
plane layer of solids from achieving densities high enough 
to trigger the gravitational instability (Weidenschilling & 
Cuzzi 1993). Even in the absence of self-sustained turbu- 
lence such as the one generated by the magneto-rotational 
instability (MRI; Balbus & Hawley 1991; Balbus & Hawley 
1998), the solids themselves can generate turbulence due 
to the backreaction of the drag force onto the gas. Such tur- 
bulence can be brought about by Kelvin-Helmholtz insta- 
bilities due to the vertical shear present in the sedimented 
layer of solids (Weidenschilling 1980; Weidenschilling & 
Cuzzi 1993; Sekiya 1998; Johansen et al. 2006), or by 
streaming instabilities induced by the radial migration 
of solids particles (Youdin & Goodman 2005; Johansen 
et al. 2006; Paardekooper 2006; Youdin & Johansen 2007; 
Johansen & Youdin 2007). In the turbulent motion, the 
solids are stirred up by the gas, forming a vertically ex- 
tended layer where the stellar gravity is balanced by tur- 
bulent diffusion (Dubrulle et al. 1995; Garaud & Lin 2004). 

But if turbulence precludes direct gravitational col- 
lapse through sedimentation, it was also shown that it al- 
lows for it in an indirect way. As solid particles concentrate 
in high pressure regions (Haghighipour & Boss 2003), the 
solids-to-gas ratio can be enhanced in the transient tur- 
bulent gas pressure maxima, potentially reaching values 
high enough to achieve gravitational collapse. Numerical 
calculations by Johansen et al. (2007) show that this is in- 
deed the case, with the particles trapped in the pressure 
maxima generated by the MRI collapsing into dwarf plan- 
ets when the gravitational interaction between particles is 
considered. They also show that the MRI is not necessarily 
needed, since the weak turbulence brought about by the 
streaming instability itself can lead to enough clumping 
under certain conditions. Another way of achieving high 
enough densities for gravitational collapse of the solid 
layer was shown by Rice et al. (2004) and Rice et al. (2006), 
where meter-sized solids concentrate prodigiously in the 
spiral arms formed in marginally gravitationally unstable 
circumstellar disks. 

Such models, however, ignored the possibility of frag- 
mentation of particles upon collisions. As the turbulence 
enhances the velocity dispersion of solids, destructive col- 
lisions become more likely. Moreover, upon destruction, 
the smaller fragments are tightly coupled to the gas and 
therefore dragged away from the midplane (Johansen et 
al. 2008), reducing the effective amount of solid material 
available for collapse. Such problem is particularly severe 
in the high mass disks investigated by Rice et al. (2004) and 
Rice et al. (2006), where the typical speeds of the boulders 
upon encounters are comparable to the sound speed. 

The fragmentation problem could be avoided if the 
accumulation of solids happened, for instance, within a 
protective environment where the collisional speeds are 
brought down to gentler values. Anticyclonic vortices 
(Marcus 1990) have been shown to favor planet forma- 
tion (Barge & Sommeria 1995; Tanga et al. 1996; Bracco et 
al. 1999; Chavanis 2000; Johansen et al. 2004) since, by ro- 



tating clockwise in the global counterclockwise Keplerian 
flow, they enhance the local shear and induce a net force on 
solid particles towards their centers. Klahr & Bodenheimer 
(2006) further argue that anticyclonic vortices would be 
less turbulent than the ambient gas, which in turn would 
lead to velocity dispersions that are low enough to prevent 
fragmentation. Vortices in disks can be the result of the 
baroclinic instability (Klahr & Bodenheimer 2003; Klahr 
2004; Petersen et al. 2007), the Rossby wave instability 
(Lovelace et al. 1999; Li et al. 2000; Li et al. 2001) or, per- 
haps, the MRI (Fromang and Nelson 2005). 

In this paper, we focus on vortices generated by the 
Rossby wave instability (RWI), which is a global instabil- 
ity where azimuthal modes experience growth in the pres- 
ence of local extrema of a quantity interpreted as a com- 
bination of entropy and potential vorticity. In the linear 
phase, the instability launches inertial-acoustic waves. The 
non-linear saturation is achieved when the Rossby waves 
break and coalesce into anticyclonic vortices. It was shown 
by Varniere & Tagger (2006) that a favorable profile of 
the entropy-modified vorticity naturally arises if the disk 
has a slow-accretion zone, such as in the layered accretion 
model of Gammie (1996). In this model, ionization is pro- 
vided by collisions in the hot inner regions, and by cos- 
mic rays in the outer disk where the column densities are 
low (a standard value for the penetration depth of cosmic 
rays is a gas column density of 100 gem -2 ). Throughout 
most of the midplane, however, the temperatures are too 
cold and the column densities are too thick for ionization 
to occur either way. The result is that, when threaded by a 
weak magnetic field, the disk displays MRI-active regions 
in the ionized layers, and a MRI-dead zone in the neutral 
parts around the midplane (Gammie 1996, Miller & Stone 
2000; Oishi et al. 2007). Matter flows towards the star due 
to the high turbulent viscosity of the MRI-active layers, 
but upon hitting the border of the dead zone, it reaches 
a region of slow accretion and the flow stalls. However, as 
the flow proceeds unabridgedly from the outer active re- 
gions, a surface density maximum forms, which triggers 
the growth of the RWI. 

The implications of this scenario for planet forma- 
tion were first explored by Inaba & Barge (2006), who 
use the RWI-unstable dead zone model of Varniere & 
Tagger (2006) to study the accumulation of solids inside 
the Rossby vortices. They confirm that the vortices are effi- 
cient particle traps, since the solids-to-gas ratio was raised 
by at least one order of magnitude, modeling the solid 
phase of the disk as a fluid. Such approximation requires 
that the size of the solid particles be much smaller than 
the gas mean free path. Since in the Minimum Mass Solar 
Nebula (MMSN; Weidenschilling 1977b) at 5.2AU the par- 
ticles subject to maximum drift have a size comparable to 
the mean free path, the sizes that a fluid approach can 
handle correspond to too strong friction, thus ultimately 
underestimating the trapping performance of the vorti- 
cal motion. In Lyra et al. (2008b; hereafter LJKP08), we 
took the works of Varniere & Tagger (2006) and Inaba & 
Barge (2006) one step further by including gravitation- 
ally interacting centimeter and meter size solids treated as 
Lagrangian particles. In that Letter, we showed that the 
solids concentrated in the vortices triggered by the RWI 
rapidly reach critical densities and undergo collapse into 
rocky planets. The resulting burst lead to the formation of 
20 rocky protoplanetary embryos in the mass range 0.1-0.6 
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Me, along with hundreds of smaller bodies following a 
mass spectrum of power law — 2.3±0.2. 

In this paper we further detail the method used in 
LJPK08, also presenting a number of new results. In the 
following section we present the dynamical equations, fol- 
lowed by an in-depth analysis of the vortices in Sect. 3. In 
Sect. 4 we analyze the formation and evolution of the pro- 
toplanetary embryos, focusing on stability against erosion 
(Paraskov et al. 2006; Cuzzi et al. 2008) and tides from the 
gas, which we identify as an important disrupting agent. 
In Sect. 6 we investigate the response of the RWI to effects 
not considered in the original analysis of Lovelace et al. 
(1999) and Li et al. (2000), and in Sect. 7 we present a dis- 
cussion of the limitations of the model. A summary and 
conclusions are presented in Sect. 8. 



2. The model 

2.1. Dynamical Equations 

We work in the thin disk approximation, using the verti- 
cally integrated equations of hydrodynamics 



dt 

du 
dt 

dx p 

~dt 
dvp 

~dt 
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where C^ ps and C^ tk are the coefficients of Epstein and 
Stokes drag, respectively. These are 
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where Ma = |Az?|/c s is the Mach number, Re = 
2a.p g \Av\/fi is the Reynolds number of the flow past the 

particle, and ]i = yf&] rcp g c s \/3 is the kinematic viscosity 
of the gas. We defer the reader to Lyra et al. (2008c) for fur- 
ther details of the drag force. The self -gravity solver is also 
explained in that paper. 

2.2. Initial and boundary conditions 

In this paper, we use the Pencil Code 1 in Cartesian and 
cylindrical geometry. The cylindrical runs do not include 
the gravity of the particles and were therefore only used 
for tests or runs without particles, as will become clear in 
the next sections. A Cartesian box was used for the pro- 
duction runs. The Cartesian box ranges x,y £ [—2.0,2.0]^. 
The resolution is 256x256, unless stated otherwise. The 
cylindrical grid ranges r £ [0.3,2.0]ro, with 2rc coverage 
in azimuth. 

The density profile follows the power law Z g =ZQr~ 03 
and the sound speed is also set as a power law c s = 



c s r 



-0.5 



. These slopes are chosen to cancel a radial de- 
pendency on the mass accretion rate m = 3nvZ g if the 
viscosity follows the Shakura-Sunyaev parametrization, 
v=#c 2 I O. The velocity field is set by the condition of cen- 
trifugal equilibrium 



r 



nl + - 

r 



1 dP d0 sg 
Za dr dr 



(12) 



We use units such that Tq=Zq=GMq=1. We choose 
In the above equations G is the gravitational constant, c SQ = 0.05 and ro=5.2 AU. We impose that the disk should 



^ g and Zp are the vertically integrated gas density and 
bulk density of solids, respectively; u stands for the veloc- 
ity of the gas parcels; Xp is the position and Vp is the veloc- 
ity of the solid particles, P is the vertically integrated pres- 
sure, c s is the sound speed, the gravitational potential, v 
the viscosity, and S the rate-of-strain tensor. The functions 
fu{Z g ) and f v (u,Z g ) are sixth order hyperdiffusion and 
hyperviscosity terms that provide extra dissipation near 
the grid scale, explained in Lyra et al. (2008a). They are 
needed because the high order scheme of the Pencil Code 
has too little overall numerical dissipation. 

The function fa is the drag force by which gas and 
solids interact. In Eq. (8), p. is the internal density of a solid 
particle, a. its radius, and Av = v p — u its velocity rela- 
tive to the gas. Cp is a dimensionless coefficient that de- 
fines the strength of the drag force. We use the formula of 
Woitke & Helling (2003) that interpolates between Epstein 
and Stokes drag 



have 30 Me of solid material within the modeled range. 
We also use a surface density of solids following the same 
slope of the surface density Ep=eEQr~ 03 , where e is the ini- 
tial solids-to-gas ratio. For the interstellar ratio e=10 -2 , the 
surface density at rg is Zq =277 gem -2 (which we round to 
Zq=300 gem -2 ). The gas disk thus has a mass of 1O _2 M 
within the modeled range. We further discuss the disk 
mass in Sect. 7.5. The gas has a Toomre Q parameter of 
30 at rg, and is throughout stable against gravitational in- 
stability. 

The dead zone is modeled as static viscosity jumps fol- 
lowing arc-tangent profiles 



tanh 



n 



Ar 



tanh 



Ar 



(13) 



1 See http://www.nordita.org/software/pencil-code 
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where r\=0.6 and r2=1.2 are the locations of the jumps 
and Ar its width. We adopt Ar=10 -2 , which provides a 
smooth jump over two grid cells only. The jump is thus 
close to a Heaviside function yet still differentiable. We 
use Vo=2.5 x 10 -5 in code units, corresponding to an alpha 
value (Shakura & Sunyaev 1973) of oc=v O / c 2 ^10 -2 . 

For the solids, we use 10 5 or 4 x 10 5 Lagrangian nu- 
merical particles. For a gas mass of 10 _2 M© and the inter- 
stellar solids-to-gas ratio of 10 -2 , each numerical particle 
therefore is a super-particle containing (in the lower reso- 
lution case) 10~ 9 M© ~ 2.7 x 10~ 2 M Moon of material. We 
use particles of radii a.=l, 10, 30, and 100 cm, as also used 
in Lyra et al. (2008c). For our nebula parameters, maxi- 
mum drift occurs for particles of 30 cm, as detailed in that 
paper. 

The particles are initialized as to yield a surface density 
following the same power law as the gas density, and their 
velocities are initialized to the Keplerian value. We use 
reflective boundaries for the cylindrical grid and frozen 
boundaries for the Cartesian grid (Lyra et al. 2008a). Both 
use the buffer zone described in de Val-Borro et al. (2006) 
to damp waves before they reach the boundary. Particles 
are removed from the simulation if they cross the inner 
boundary. 



3. Vortices 

The trapping mechanism of vortices is not only due to its 
being a high-pressure region, but mainly due to the vor- 
ticity of the flow. In an anticyclonic vortex (cyclonic vor- 
tices are destroyed by the Keplerian shear), the motion oc- 
curs in the same sense as the local shear, i.e., the gas ro- 
tates clockwise. Therefore, at the antistellar point the angu- 
lar momentum is decreased with respect to a non-vortical 
flow; and conversely increased at the substellar point. As 
a result, the gas at the antistellar point is accelerated in- 
wards, while the gas at the substellar point is accelerated 
outwards. A net centripetal force towards the eye ensues. 
The streamlines of vortices (or vortex lines) are a set of 
Keplerian ellipses with the same semimajor axis but dif- 
ferent eccentricities, being circular in the center and more 
eccentric outwards (Barge and Sommeria 1995). We show 
contours of | u | on the surroundings of one of the giant vor- 
tices, in the upper left panel of Fig. 1. As the gas drags the 
particles, the particles also revolve around the vortex eye. 
But because the gas-solids coupling is not perfect, the par- 
ticles lose angular momentum and sink deeply towards 
the center. In the next subsections we describe some of the 
properties of the vortices present in our simulations 




Fig. 1. Enlargement around one of the vortices in a snapshot at 
20 orbits. Contours of \u\ are superimposed on the gas surface 
density plot, showing that the density enhancement is associated 
with intense vorticity. In the upper panel we show the multi- 
phase (total) surface density of solids, whereas in the middle and 
lower panels we show the contribution of each particle species. 
The vortical motion preferentially traps particles of # =1O and 30 
cm. 



The quantity T is defined as 

ZD 

= k 2 -Acv 2 -c 2 / (L s L p ) 

where 



3.1. Launching Mechanism - the RWl 

The vortices in LJKP08 are triggered by the Rossby wave 
instability (RWI), a case of purely hydrodynamical insta- 
bility in accretion disks. Considering azimuthal perturba- 
tions to the inviscid Euler equations, Lovelace et al. (1999) 
and Li et al. (2000) find that instabilities exist when the fol- 
lowing quantity has a local extremum 
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Fig. 2. Evolution of the pressure bump with different widths of 
the viscosity jump (uppermost panels; the dashed line is the 
Heaviside jump, for comparison). The violation of the Solberg- 
Hoiland criterion (lowermost panels) is a conservative indica- 
tion that the threshold of instability of the RWI was reached. As 
the width of the jump increases, the threshold takes increasingly 
longer time to be breached. For jumps smoother than two scale 
lengths, the threshold was not reached (up to 500 orbits). 



are the radial length scale of the entropy and density vari- 
ations, respectively. 7 is the adiabatic index. For corota- 
tional modes (Aco = to — mQ <C k) in a barotropic (L s — » 
00) disk, the quantity T reduces to ZOk~ 2 , which is read- 
ily identified with (half) the inverse of vortensity £ 



co z /L 
| V x u\ z 

r dr ( r ^) 



(19) 
(20) 



K 

2fi' 



which in turn led Lovelace et al. (1999) to interpret C as an 
entropy-modified version of, or generalized, potential vor- 
ticity. An extremum in the profile of C can be generated, 
for example, by a pressure bump somewhere in the disk. 
The dispersion relation of the disturbances is analogous to 
the dispersion relation of Rossby waves in planetary atmo- 
spheres, hence the name of the instability. 



3.2. How sharp need the jump be? 

Although formally the instability is triggered by a min- 
imum or maximum of C, ones finds in practice that the 
amplitude and radial width of the pressure bump present 
critical values beyond which instability does not occur. Li 
et al. (2000) find that typically, a pressure variation of 10- 
20% over a length similar to the scale height of the disk is 
sufficient to trigger the instability. The threshold, however, 
is problem dependent, depending - among other things 
- on the geometry of the pressure variation (step jump 
or Gaussian, for example; see Li et al. (2000) for details). 
Due to this, we use a more general criterion to assess the 
threshold of instability. Li et al. (2000) note that the thresh- 
old of instability for the RWI is always reached before the 
Solberg-Hoiland criterion for stability of axis-symmetric 
disturbances 

k 2 + N 2 > (21) 

is violated. Eq. (21), therefore, provides a conservative es- 
timate of whether or not the RWI is excited. The Solberg- 
Hoiland criterion is easily understood. In a pressureless 
disk, the condition k 2 > suffices to determine stability. 
The other term 



N z 



1 dP 



1 dZ 
Z~dr 



1 dP\ 
jP~dr) 



(22) 



is the square of the Brunt-Vaisala frequency, associated 
with the oscillations of buoyant structures in the pres- 
ence of an entropy gradient. Physically, Eq. (21) means 
that a mode that is unstable /stable to shear can be stabi- 
lized/ destabilized by pressure gradients and vice-versa. 

We measure the epicyclic and Brunt-VaisaTa frequen- 
cies in a series of ID simulations where we varied the 
width Ar of the viscosity jump (Eq. (13)). We find that for 
locally isothermal simulations, the Solberg-Hoiland crite- 
rion is violated at the outer edge for Ar < 0.04, which 
is slightly less than one scale height. In these isother- 
mal simulations, the criterion depends almost solely on 
the epicyclic frequency because, as the temperature does 
not rise with compression, the pressure does not change 
enough for N 2 to go appreciably negative. 

To assess the effect of non-barotropic behavior, we re- 
place the locally isothermal flow by an isentropic one with 
an adiabatic equation of state 



as 

dt 



-(u-V)S+f x {S) 
Ec 2 / 7 



(23) 
(24) 



which means that we allow heating and cooling by com- 
pression and rarefaction only, excluding viscous heating 
and radiative cooling. In Eq. (23), S = In P — 7 In Z is the 
vertically integrated specific entropy of the gas. The func- 
tion f x (S) is a sixth order hyperconductivity term, analo- 
gous to hyperdiffusion for density. The adiabatic index is 
7=7/5. 

The results are illustrated in Fig. 2, where we plot the 
viscosity profile (upper panel) for different widths Ar, and 
the time evolution of the density, pressure and temper- 
ature bumps. The lower panels measure if the Solberg- 
Hoiland criterion was violated. We find that the pressure 
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bump sharpens considerably compared to the isothermal 
case, due to the high temperatures associated with the 
compression. The consequence is that the Solberg-Hoiland 
criterion is violated by viscosity jumps up to Ar < 0.12, 
i.e., 3 times broader than in the isothermal simulations. In 
this non-isothermal case, it is mostly the Brunt- Vaisala fre- 
quency that leads k 2 + N 2 to negative values. The effect of 
increasing the width is mainly of slowing the evolution of 
the quantity k 2 + N 2 towards negative values. It only takes 
five orbits for Ar = 0.01, but it takes 350 orbits when Ar is 
increased to 0.1. In Fig. 2 we state Ar in terms of the scale 
height H=0.05ro- We present a resolution study of vortex 
excitation in Appendix A. We also address the issue of vor- 
tex survival in a non-static dead zone in Appendix B. 

3.3. Steady-state dead zone 

If no transport happens in the dead zone, matter can do lit- 
tle more than piling up there as the inflow proceeds from 
the active layers. However, the accumulation of matter 
cannot proceed indefinitely since, as matter piles up, the 
conditions for gravitational instability would eventually 
be met (Armitage et al. 2001). The gravitational turbulence 
that ensues (Lodato 2008) would therefore empty the dead 
zone as the excess matter accretes, thus re-starting the cy- 
cle. 

However, local simulations show that the dead zone 
has some level of residual turbulence. This happens be- 
cause the turbulence on the active layers induce small lev- 
els of Reynolds stress in the dead zone (Fleming & Stone 
2003). If the inertia of the midplane layer is not too high 
(Oishi et al. 2007), this forced turbulence can lead to mod- 
erate oc values with non-negligible transport 2 . 

Terquem (2008) shows that steady state solutions in ID 
models exist in this case, as the dead zone gets denser and 
hotter to match the condition of constancy of the mass ac- 
cretion rate with radius, d r (vE)=0. In this case, the steady 
state will have an vj viscosity value in the active layers 
and a lower in the dead zone. 

Vortex formation by the RWI requires the presence of 
a pressure maximum. In our model, and that of Varniere 
& Tagger (2006) and Inaba & Barge (2006), the pressure 
maximum comes about by stalling the accretion flow in the 
border of the dead zone. There is no requirement that the 
dead zone should have zero viscosity, just a viscosity sig- 
nificantly lower than that of the active regions. We tested 
different values of / vj, and found that changing it from 
to 0.1 has little effect on the instability. For higher values, 
the Solberg-Hoiland criterion takes increasingly longer to 
be violated. For Vd/vj = 0.5, the Solberg-Hoiland cri- 
terion is violated after 60 orbits. We also notice that the 
steady-state dead zones of Terquem (2008; see Fig. 3 of that 
paper) have the surface density varying by more a factor 
of ^10 over a few scale lengths at the inner edge. Such 
profiles violate the Solberg-Hoiland criterion, so the RWI 
is expected to be excited in those conditions as well. 



2 Another alternative is local ionization provided by the de- 
cay of the short-lived radioactive nuclide 26 Al (Umebayashi & 
Nakano 2009), although Turner & Sano (2008) show that the free 
electrons given out by this low ionization source would quickly 
recombine on the surface of ^m-sized dust grains. 



With drag force Without drag force 




5 10 15 20 25 30 5 10 15 20 25 30 

t 




5 10 15 20 25 30 5 10 15 20 25 30 

t 



Fig. 3. ID simulation of collapse of 1000 particles. With drag force 
the kinetic energy of the particles is efficiently dissipated and the 
particles collapse at subgrid scale towards infinite density. When 
the drag force is excluded the system cannot dissipate energy and 
a halo of particles, 10 grid cells wide, is formed. 



0.8 




Fig. 4. Evolution of the most massive clump formed (solid line), 
which we traced back in time from the end of the simulation. 
It differs from the instantaneous most massive clump (dashed 
line) because the clumps have different feeding rate and can also 
experience mass loss, as in the episode that happened at w 90 
orbits (see text). 



4. Embryos 

4.1. Drag force cooling and compactness 

The embryos formed in our simulations present a number 
of interesting features. We first would like to address the 
issue regarding their physical size. The embryos consist 
of a cluster of a large number of particles, held together 
by their collective gravitational pull. But are they strongly 
bound like solid objects? Or do they consist of loosely cou- 
pled objects in the same potential well? To answer this 
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question we measure the rms spatial dispersion of the par- 
ticles inside the cluster, defined as 



El 

i=l 



'CM I 



(25) 



where n is the number of particles within the Hill sphere of 
the clump, r is the vector radius of each particle and r CM is 
the vector radius of their center of mass. We take this value 
as a measurement of the "radius" of the embryo. The most 
massive embryo has a radius of 1.13 R@. This compactness 
corresponds to a tenth of a thousandth of the resolution 
element. 

Such compactness is due to the dynamical cooling pro- 
vided by the drag force. We illustrate this in Fig. 3. The fig- 
ure shows the results of ID simulation with a thousand in- 
teracting particles with and without gas drag. Without gas 
drag the particles have no means to dissipate energy and 
perform oscillations about the center of mass. The very in- 
ner particles show virialization, while the outer particles 
form a halo extending for a radius of 10 grid cells in aver- 
age. 

When including gas drag, the system gets so dissipa- 
tive that the kinetic energy is soon lost and the ensemble of 
particles collapses. The exponential decay of the particles' 
rms position seen in the upper left panel of Fig. 3 shows no 
sign of flattening, down to a millionth of a resolution ele- 
ment. This leads us to infer that collapse to zero volume is 
ongoing. This is of course expected, since no mechanism to 
provide support against the gravitational pull is present. 

In view of this, the question is why our planets, which 
are subject to drag forces, do not shrink to zero as well but 
stabilize at a very small but finite radius. We are drawn to 
the conclusion that this is a numerical issue. The tests of 
Fig. 3 were done with a fixed time-step. But when the par- 
ticles cluster together to form a planet in our simulations, 
they end up dominating the time-step. The position up- 
date x(t) = xq + vAt therefore occurs with the maximum 
At allowed by the Courant condition, which is that the 
fastest particle should move by one grid cell. Due to this, 
the time resolution of the subgrid motion around the cen- 
ter of mass of the cluster is under-resolved. With the over- 
shot At, the particles that are attracted towards the center 
of mass of the clump will end up in a position past it. In the 
next time step they will be attracted to the center of mass 
from the other side, but will once again overshoot it. The 
result of this is that the particles will execute undamped 
oscillations, leading to a finite rms radius. We performed 
tests like those of Fig. 3 with a particle-controlled variable 
timestep, confirming this explanation. We conclude that 
the fact that the most massive embryo has a stable rms 
radius compatible with its mass is but a deceptive coin- 
cidence. 

We stress that this drag force cooling will cease to be ef- 
ficient as the solids-to-gas ratio grows too large (pp/pg ^> 
1), because in this case the backreaction would be too 
strong and the gas would simply be dragged along with 
the particles. In this case, a Keplerian disk of solids might 
form, accreting matter onto the planet due to collisions be- 
tween the orbiting solids. This accretion regime is never- 
theless beyond the scope of the current paper. 
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Fig. 5. Time series of the mass loss episode. Due to gravitational 
tides from the massive vortex whence the embryo formed, a large 
chunk of particles was detached from the original cluster. At 87 
orbits, a separation 4 times the Earth-Moon distance is seen. The 
separation grows and two orbits later the two bodies do not over- 
lap Hill radii, thus counting as separate objects. 



4.2. Mass loss 

In Fig. 2a of LJKP08 we showed the evolution of the most 
massive clump. However, as the clumps have different 
feeding rate and some of them experience mass loss, the 
most massive clump shown there is not always the same 
clump. In Fig. 4 we contrast this with the evolution of the 
most massive clump at the end of the simulation, which we 
tracked backwards in time. Such clump started in the inner 
disk, showing 0.8 Mj^ars by 40 orbits. By this time, the most 
massive clump was a 3 MMars clump in the outer disk. 

The most remarkable feature of this plot is the mass 
loss event at 90 orbits. Fig. 5 shows that it consists of the 
detachment of a 0.8 M^iars object from the original cluster, 
of 6.7MMars- The detachment is already seen at 87 orbits, 
although the separation is quite small (4 times the Earth- 
Moon mean separation, -Rem)- At 89 orbits, the smaller ob- 
ject left the Hill sphere of the more massive embryo. They 
finally appeared as separate objects, and the maximum 
mass decreased. 
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Fig. 6. The state of the disk before (a-d) and after (e-h) the mass loss episode. The conspicuous difference between them is that the 
embryo has left its parental vortex from one snapshot to the other. It is seen as a bright spot in panels / and h, at (X,Y)=(-0.65,-0.19). 
In panel b (before the mass loss), the embryo is at (X,Y) =(-0.40,-0.53) but not easily spotted among the swarm of solids inside the 
vortex. Panel i shows a horizontal slice through this location, in which we see that the density of solids does not peak much higher 
than the gas density at the location of the embryo (panel ;). Significant gas tides are expected as the gravitational potential (panel k) 
and acceleration (panel I) have similar contributions from the gas and solid components. 



We see evidence that this puzzling behavior is due to 
gravitational tides from the gas. The gas is too pressure- 
supported to undergo collapse, but the vortices concen- 
trate enough material to yield a considerable gravitational 
pull. We illustrate this in Fig. 6, where we show the state of 
the disk before the mass loss episode (at 80 orbits, Fig. 6a- 
Fig. 6d) and after that (at 100 orbits, Fig. 6e-Fig. 6h). The 
plots show the surface densities of gas and solids, and the 
potential associated with them. Even though the clump- 
ing of solids yield a considerable gravitational pull (Fig. 6d 
and Fig. 6h), these figures show that the dominant contri- 
bution to the gravitational potential comes from the gas - 
more specifically from the vortices, where the gas density 
peaks one order of magnitude denser than the initial con- 
dition. 

The most massive clump is located in the inner disk at 
(X,Y)=(-0.40,-0.53) in Fig. 6b, not clearly identifiable amidst 



the other particles trapped inside the vortex. However, 
the embryo is immediately observable as the bright point 
at (X,Y)=(-0.65,-0.19) in Fig. 6h (also visible in Fig. 6f, al- 
beit less prominently). There are two features in this plot 
that are worth noting. First, by comparing the location of 
the embryo in these figures with the location of the vor- 
tices, we notice that the planet has left its parental vortex. 
Second, the inner vortices have undergone the transition 
from the m=3 to the m=2 mode. Due to merging, their gas 
density has increased, with dramatic consequences for the 
embryos within them. 

We assess how the gravity of the gas influences the mo- 
tion of the particles (Fig. 6i-Fig. 61). In Fig. 6i we take a 
horizontal density slice at the position of the most mas- 
sive embryo at 80 orbits. Fig. 6j is an enlargement of Fig. 6i 
around X=-0.53, where the embryo is located. We see that 
the densities of solids and gas peak at similar values. The 
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next figures show the gravitational potential (Fig. 6k) and 
acceleration (Fig. 61) around the embryo. The gas produces 
a deeper gravitational well, albeit smoother than the one 
displayed by the solids. In the acceleration plot it is seen 
that the pull of the gas is greater than the pull of the em- 
bryo already at a distance of just 0.26 AU (0.03 in code 
units, corresponding to two grid cells) away from the cen- 
ter. And even where the pull of the solids is strongest (one 
grid cell away from the center), the gravity of the gas still 
is an appreciable fraction of the gravity of the solids. Tides 
from the gas are unavoidable. 

It is beyond the scope of this paper to consider the full 
mathematical details of the theory of tides, especially be- 
cause the two bodies (the vortex and the embryo) are ex- 
tended. Instead, we consider the following toy model. The 
tidal force F T experienced by the planet is proportional to 
the gradient of the acceleration a induced by the gas. It is 
also proportional to the radius R of the planet: Fj oc R Va. 
Since Va = — V 2 <£, according to the Poisson equation, the 
tidal force is proportional to the local value of the density 
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F T oc Rp g . 



(26) 



We consider the 3D volume density to avoid the 
requirement of using the Dirac delta in the 2D case. 
Considering the planet spherical, Newton's second the- 
orem holds and we can write Fq = — GM/R 2 
for the planet's (self-) gravitational force at its surface. 
Substituting M = 4/3nppR 3 , we have Fq oc Rp v , so 



r F T Pg 
b = — a / 
F G P V 



(27) 



i.e., the ratio of the disrupting tidal stresses to the self- 
gravitating forces that attempt to keep the planet together 
is directly proportional to the gas-to-solids ratio. At 80 or- 
bits, as seen in Fig. 6j this ratio is around unity. As the 
vortices undergo the transition from the m=3 to the m=2 
mode, their peak density increases (while the planet re- 
mains at a constant mass), the tides eventually become 
strong enough to cause the mass loss event of Fig. 5. This 
effect will probably be less dramatic in 3D simulations be- 
cause, as the particles settle in the midplane, the ratio of 
the volume gas density to the bulk density of solids Pg/pp 
is expected to be much lower than the ratio of the column 
gas density to the vertically integrated surface density of 
solids Zg/Zp. 

4.3. Erosion? 

Cuzzi et al. (2008) points that erosion is of prominent im- 
portance in the stability of self-gravitating clumps of par- 
ticles. They put forth a model where self-gravity plays a 
role analogous to that of surface tension in liquid drops, 
preventing disruption against ram pressure forces from 
the gas. The clumps are held together by self-gravity if the 
gravitational Weber number (in analogy with the surface 
tension case) is less than a critical value, close to 1. The 
gravitational Weber number is defined as the ratio of the 
drag to self-gravitational accelerations 
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Fig. 7. Can self-gravity alone maintain the formed embryos to- 
gether? We switch off the hydrodynamics at the last snapshot 
and run for additional 200 orbits. The most massive embryo, 
which had already left the parental vortex (at « 90 orbits), does 
not disperse. 



Cuzzi et al. (2008) further point that in numerical mod- 
els, artificial viscosity can largely exaggerate the disrupt- 
ing effect of the ram pressure. This happens because, as the 
clumps are small, they are deeply within the viscous range 
of the grid, whereas in the real solar nebula the dissipation 
happens at much smaller scales. The Reynolds number of 
the flow past the particles is therefore much smaller than 
what a real clump would experience 3 , and the exagger- 
ated viscous stresses might de-stabilize potentially stable 
clumps. 

It is interesting to assess if this mechanism plays a sig- 
nificant role in our models, or even if it can account for at 
least some of the mass loss events. 

We can estimate how important erosion will be for our 
clumps the following way. We approximate the clumps as 
single point masses so that | V<£ S g| ~ GM/r 2 , where M is 
the total mass of the clump. Plugging this in Eq. (30), we 
write the gravitational Weber number as 



We G 



3pC D \Av\ 2 r 2 
SGMa.p. 



(29) 



For Epstein drag (Eq. (10)), Co does not depend on a m . 
So, for all other quantities being constant, we expect Wee 



3 The Reynolds number of the flow past a clump can be written 
as Re = Rrms^rms/v. At the grid scale our choice of viscosity is 
usually 3 x 10 17 cm s -1 (it decreases very fast as we go to smaller 
wavenumbers, as k 6 ). For a clump as the ones of this study, of 
R rms =10 4 km and v rms =l ms _1 , the Reynolds number is Re « 
3 x 10~ 7 . At such incredibly low Reynolds numbers, inertia plays 
no role. The self-gravity of the particles, therefore, is not holding 
the cluster together against drag forces from the gas, but against 
largely exaggerated viscous stresses. In comparison, in the solar 
nebula, the molecular viscosity is much lower and the Reynolds 
number is expected to be > 10 6 (Cuzzi et al. 2008). 
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to decrease linearly as a. increases. In other words, self- 
gravitating clumps of larger particles should be more sta- 
ble than clumps composed of smaller particles. 

We can simplify the Wee by writing the mass M as 
M=Tcr^ ms E pr and the drag force as |/^|=Mac s /T. Thus, at 



Relative vs Collisional Velocities 



r=r r 



We G = 



Ma c s 



(30) 



We confirm in Appendix C that Eq. (30) is sufficiently 
accurate in predicting the onset of erosion. For our choice 
of parameters, 



We G 



17 Ma 



(31) 



so for a flow of Ma « 10 -2 , a clump of Z p =l at r « 1 will 
be stable if rO > 0.1. The embryos we consider are formed 
predominantly of particles of 10 cm (tH^O.1) and 30 cm 
(tH^O.3). We conclude that erosion might play a role in 
our case. The sharp dependence of Wee on the distance in 
Eq. (31) also means that embryos at the inner edge of the 
dead zone are more prone to erosion than the ones at the 
outer edge. We will develop this further in Sect 5.2.2. 



4.4. Can the embryos be held together indefinitely? 

To answer the question of how long lived these clumps 
are, we take the last snapshot of the simulation and switch 
off the hydrodynamics. The particles now move under the 
influence of the stellar gravity and their own self-gravity 
only. We run for additional 200 orbits to assess is self- 
gravity alone can maintain the cluster of particles together. 
The result is shown in Fig. 7. 

The clustered particles do not disperse, and the most 
massive embryo maintains the same mass for the addi- 
tional simulating time. We do not see difference in the de- 
gree of compactness of the most massive embryo. As the 
vortices are shut down, the unbound 1 cm sized particles 
that were too well coupled to the gas to be dragged into 
the eye are released from they vortical confinement and 
spread over a wider annulus. 

We have no reason to suspect that the situation will 
change over longer timescales. We conclude that the em- 
bryos can be held together indefinitely. 

4.5. Collisions 

As stated before, one of the problems that solids accumu- 
lation inside vortices can potentially solve is the issue of 
fragmentation of particles upon collisions, a drawback for 
both coagulation (Brauer et al. 2008a) and gravitational in- 
stability models (Rice et al. 2006, Johansen et al. 2007) of 
planetesimal growth. 

To assess if fragmentation poses a significant barrier for 
the formation of the protoplanetary embryos in this study, 
we take the most massive embryo by the end of the sim- 
ulation, flag the particles that constitute it and trace them 
back in time, calculating their collisional velocity history. 
The collisional speed for each particle is calculated by tak- 
ing the closest neighbor to that particle within the range 
of a grid cell. There is, however, a subtlety concerning the 
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Fig. 8. The relative velocities between particles are contrasted 
with collisional velocities calculated from Eq. (32). At the initial 
condition (t=0) the relative velocities between particles follow 
what is expected from the Keplerian shear (dashed line), reach- 
ing as far as 300 m s -1 inside a grid cell (which is 31 Rem wide). In 
contrast, the true collisional velocities are much lower. The figure 
shows that estimating fragmentation with the relative velocities 
would greatly overestimate its likelihood. Time is quoted in or- 
bits (t K =2n/n ). 



difference between collisional velocities and relative veloc- 
ities. A collision between particles i and ; only happens 
when the separation r z y between them tends to zero. For 
our resolution and choice of rg, a grid cell is 0.08 AU wide, 
thus existing plenty of room for subgrid Lagrangian dy- 
namics. In particular, the velocity difference due to the 
Keplerian shear between the inner and outer radial bor- 
ders of a grid cell can be significant. At the inner edge of 
the dead zone of the model presented in this paper (3.12 
AU) for instance, this difference amounts to 434 ms" 1 . As 
this velocity difference is due to the separation between 
particles, it vanishes when r z y tends to zero, thus never con- 
tributing to the true collisional speed. In Fig. 8 we show 
these uncorrected relative velocities in the initial condi- 
tion and in a snapshot at 40 orbits, plotted against sepa- 
ration. The clear correlation follows what is expected from 
the Keplerian shear (dashed line) in the initial condition. 

The gas motion adds another velocity that has to be 
taken into account. Solid particles are dragged by the gas 
motion, yet gas streamlines never intersect. The gas mo- 
tion itself thus introduce velocities that never participate 
in collisions. We correct for these two by the following pro- 
cedure. For each particle in the pair involved in a collision, 
we consider its velocity Av relative to the gas (the same 
quantity that appears in the drag force, Av=Vp — u). We 
then define the collisional velocity vector as 



"call 



Av(xt) 



Av(x 



(32) 



Equation (32) ensures that tracer particles, which fol- 
low the gas streamlines, should never experience collision. 
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Accretion history Maximum collisional velocity max(v coll ) vs. Separation 
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Fig. 9. a). Accretion history of the most massive embryo. Collapse happens at 40 orbits, and further accretion of particles happens 
through the next 18 orbits, after which the maximum mass is attained. This plot defines the time to that each particle is accreted. 

b) . Histograms of the maximum collisional velocities experienced by a given particle, before (dot-dashed line) and after (solid line) 
its accretion time to. The solid line histogram represents the maximum collisional speed from to until the end of the simulation. The 
dot-dashed line histogram covers a time interval of ten orbits before to- As the vortices form at ^30 orbits, the latter better represents 
v co \\ under vortex trapping. The vast majority of the particles are in the bin of 0-5 ms" 1 . The smaller window zooms into this bin. 
The black line corresponds to bins of lm s -1 whereas the gray line to bins of 0.1m s _1 . Both represent v co \\ after to- It shows that most 
of the particles never experience collisions more violent than 1 ms -1 . 

c) . Same as b). but plotting the maximum collisional speeds as function of separation between the particle and its nearest neighbor. 
Before to, three populations are seen, of low speed at short separations, high speeds at short separations, and high speeds at large 
separations. Only the second group would have experienced destructive collisions. After to, 99% of the particles belong to the first 
group. The correlation with distance for the first group is not due to the Keplerian shear (dashed line) or the residual shear present 
in Av (dotted line). 



Furthermore, as the shear dependence of Av=?jvk is much 
smaller than the shear dependence of v p =vk (where rj = 
(l/2)(H/r) 2 (31nP/31nr) = 3.75 x 1(T 3 is the pressure- 
correction factor), we can assume that it also corrects for 
this quantity. Figure 8 also shows the collisional speeds 
calculated by Eq. (32). The dependence on separation was 
greatly suppressed. 

The results of the collisional velocity history of the par- 
ticles that constitute the most massive embryo at the end 
of simulation are plotted in Fig. 9. 

Figure 9a shows a cumulative plot of the mass of the 
embryo, which defines the time to that each flagged par- 
ticle was accreted. We define the time of accretion as the 
moment when the particle approached the grid point x near 
nearest to the maximum of particle number density (also 

defined by the flagged particles) by less than d^^dxy/l 
(the grid cell diagonal) and kept 

\Xp — #near| ^ ^diag (33) 

until the end of the simulation. Although this is not as 
strict as the definition of accretion we have been using 
before (based on the Hill criterion and escape velocity), 
this simpler criterion captures what happens before col- 
lapse (i.e., before the maximum of particle number den- 
sity becomes a bound protoplanetary embryo) and serves 
well our purpose of illustrating the behavior of collisional 
speeds at close separations. The first episode of accre- 
tion takes place at 40 orbits, coinciding with the time that 
the clump of particles became bound (in accordance with 
Fig. 4). Further accretion proceeds over the next 18 orbits, 
with the maximum mass being attained at 58 orbits. No 
other particle was accreted after this time. 



Figure 9b shows histograms of the maximum col- 
lisional speed that a particle experienced before (dot- 
dashed line) and after (solid line) accretion. The latter 
refers to the maximum of v co \\ taken between to and the 
end of the simulation (200 orbits). The former refers to a 
time interval of 10 orbits before to- As the vortices in the 
inner disk just fully develop at ^30 orbits, the dot-dashed 
histogram is more representative of a situation where par- 
ticles are inserted in a disk with existent vortices. The con- 
clusion is striking: the vast majority of the particles that 
constitute the embryo never experienced a collision more 
violent than 1ms -1 . 

In Fig. 9c we show the maximum collisional veloci- 
ties of Fig. 9b as a function of the separation between a 
given particle and its closest neighbor, also before and af- 
ter accretion. The distribution before accretion is trimodal, 
with particles with low speeds (<lms -1 ) at small separa- 
tions, particles with high speeds (<20ms -1 ) at small sep- 
arations, and particles with high speeds at large separa- 
tions (>10^em)- Only the second group of particles would 
have undergone fragmentation. The first group is below 
the fragmentation velocity threshold, whereas the large 
separations of the third group imply they never experi- 
enced an encounter close enough to lead to a collision. 
After accretion, virtually all particles (99%) belong to the 
first group. 

In Fig. 10 we plot the collisional velocities versus sep- 
aration at selected snapshots instead of historical maxima. 
In these plots, we only used the particles for which Eq. (33) 
was satisfied, i.e., considering only the collisions that are 
participating on the formation process of the embryo. At 
30 orbits, a small number of particles is observed (87), 78% 
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Fig. 10. Collisional speeds as a function of separation at selected 
snapshots at 30, 35, 40, and 45 orbits {tj^=2n/ Oq). The figure only 
shows particles that are closer than the distance of a grid cell to 
the clump of maximum number density The group at low sepa- 
rations and high speeds observed in the upper plots would un- 
dergo fragmentation. Nevertheless, they represent only a minor- 
ity of particles (22% at 30 orbits, 8% at 35). At the time the cluster 
gets bound (at 40 orbits), the vast majority of particles (99%) is 
safely at low speeds. The presence of a fourth group at 45 orbits, 
at very low separations (<10 _3 Rem) an d with very low colli- 
sional speeds (<10~ 3 ms" 1 ) indicates that collapse towards zero 
volume is ongoing. 

of these showing safe collisional speeds (<10ms _1 ). 5 or- 
bits later 119 particles are within the grid cell of the form- 
ing embryo, 92% of which show gentle collisions. At the 
time the overdensity gets bound (40 orbits), it is formed 
by 639 particles, with just 7 of these showing collisional 
speeds greater than 10ms -1 . At 45 orbits, all 877 particles 
display safe speeds. The tendency seen in this time series 
towards an increasing number of encounters at low sepa- 
rations and at low collisional speeds indicates that collapse 
towards zero volume is ongoing. Indeed, at 70 orbits, we 
observe that most of the particles occupy the same point 
in space. 

5. Size distribution 

In our simulations, we considered the solid phase of the 
disk represented by particles of 1, 10, 30, and 100 cm. In 
this section we discuss a number of issues, relevant to the 
simulations, related to having a size spectrum instead of 
single-phasing. 

5.1. Aerodynamical sorting 

One of the most prominent features of the embryos formed 
in our models is that they are composed primarily of 



Fig. 11. Aerodynamical sorting for the particles trapped in the 
vortical motion. The figure is centered at the vortex shown in 
Fig. 1, at 18 orbits. The unit of length is the Earth-Moon mean 
distance, and the Y-coordinate points to the star. As particles of 
different size have different friction times, differential drag oc- 
curs inside the vortex, effectively sorting the particles spatially 
by size. 



same-sized particles. This is mostly due to aerodynamical 
sorting. As particles of different size have different fric- 
tion times, differential drag occurs inside the vortex, effec- 
tively sorting the particles spatially by size. Moreover, the 
stationary point is determined by a balance between the 
Coriolis and the drag force, in such a way that the eye of 
the vortex is the stationary point only for r s — > 0, or perfect 
coupling. In general, the stationary point is azimuthally 
shifted with respect to the eye, according to the particle 
size (Youdin 2008). 

The aerodynamical sorting inside the vortex can be 
seen in Fig. 11, which corresponds to the vortex of Fig. 1. 
As similar particles drift alike, streams of same-sized par- 
ticles are clearly seen in the vortical flow 

5.2. Differences between the inner and outer embryos 

Another interesting feature of our results is that once the 
vortices are formed, they easily trap the 10 cm and 30 cm 
particles both in the inner and in the outer edge of the 
dead zone. Yet, by the end of the simulation there seems 
to be a preference for the embryos in the inner disk to 
be composed of larger particles (30-100 cm), while in the 
outer disk, more embryos formed of the smaller particles 
(1-10 cm) are seen. We explain these two features in the fol- 
lowing sections. 
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5.2.1. Preferential sizes in different locations of the nebula 

The first feature (inner and outer vortices equally trapping 
particles of same size) follows from the fact that the gen- 
eral drag coefficient (Eq. (9)) yields a nearly flat profile for 
the radius of the particle with maximum drift (tQ=1) ver- 
sus distance. This is seen by calculating the stopping times 
t s of the different particles as a function of the dynamical 
variables 



t s = rO 

= V32^ \p. (Kn' + l) 2 
Kn'Ma Z g (Kn' 2 c|f + C* k ) 

and calculating the radius for a given stopping time 
t s . Substituting C^ ps =16/(3Ma) VT/ri (Eq. (10) at sub- 
sonic regime) and Stokes law C^ tk =24/Re (with Re = 
3Ma/Kn\/ zr/8) yields the quadratic equation for a. 

2a 2 . + 3\a. ^ = (35) 

np. 

and, as the radius is positive, the solution is unique 

3A / r 16Z g t s \ 

a - = T{f+3^;- 1 ) (36) 

Fig. 12 shows this curve as a function of radius for 
t s =l and our initial parameters. The particle radius of 
t s =l predicted by both pure Epstein and pure Stokes drag 
are shown for comparison. The mean free path A is also 
shown. The figure shows that the curve is quite flat for 
our choice of parameters, so inner and outer vortices must 
have similar efficiency on trapping particles of a given 
size. 

5.2.2. Tidal disruption and erosion of the inner embryos of 
a. =10 cm 

We see evidence that the second feature (the absence of 
clumps of 10 cm in the inner disk at later times even 
though they were formed) is due to tidal disruption in 
the same episode that lead the most massive embryo to 
lose mass. We illustrate this in Fig. 13, a time series of the 
gas and the solid phases, the latter split into the 4 differ- 
ent particle sizes. The upper plots, at 75 orbits, show that 
embryos of 30 cm and 10 cm were formed in both the inner 
and outer vortices. 

The gas plots of the time series of Fig. 13 illustrate the 
transition from the m=3 to the m=2 undergone by the inner 
vortices, as mentioned in Sect. 4.2. This raises their den- 
sity so the tides get stronger. The plots of the 30 cm phase 
show that the embryos composed of these particles split 
into smaller objects, which nevertheless can still keep their 
physical integrity. 

The fate of the clusters composed of particles of 10 cm 
is different, though. As the gas density increases, so does 
the gravitational Weber numbers of the embryos. Erosion 
starts to play a more significant role. As the density in- 
creases inside the inner vortices (a factor 5 relative to the 
initial condition at 75 orbits; 8 at 200 orbits), the embryos 



of 10 cm particles start to behave more and more like the 
rf^O.01 clusters of Fig. C.l. At 85 orbits, one of embryos 
of 10 cm particles in the inner edge was destroyed. At 95 
orbits, a second embryo was disrupted. At 105 orbits, the 
third embryo was also destroyed by the combined effect 
of tides and erosion. 25 orbits later, the 10 cm particles 
have dispersed through the inner edge of the dead zone. 
The tides from the gas prevent them from assembling once 
again. 

The outer vortices never get as strong as the inner ones. 
The result is that although the inner embryos of 10 cm par- 
ticles are destroyed, the outer ones are kept until the end 
of the simulation. 



6. The response of the RWI to the drag 
backreaction and self-gravity 

The evolution of the RWI was studied analytically for the 
case of a low mass dustless disk only. In Fig. 14 we show 
how the effects of gas self-gravity and backreaction from 
the solids affect the evolution of the instability. 

The upper panels of Fig. 14 show a disk without solids 
and without self-gravity. In the middle panels we included 
self-gravity, while in the lower panels we included solids. 
The appearance of the disk in the three simulations is 
shown in selected snapshots at 5, 10, and 15 orbits. 

The self-gravitating and non-selfgravitating dustless 
cases (upper and middle panels) look similar, with the 
RWI being excited first in the outer edge of the dead zone. 
However, there is a crucial difference between them. The 
snapshot at 15 orbits shows a prominent m=2 mode in 
the outer disk for the non-selfgravitating case, while the 
run with self -gravity displays a dominant m=5 mode at 
the same time. This puzzling result is made even more in- 
teresting by recalling that the gas is gravitationally stable 
(Q^30). That the growth rates of different modes vary that 
significantly for such a value of Q is indicative that the dis- 
persion relation of the RWI is probably remarkably sensi- 
tive to self-gravity. 

The simulation with drag backreaction (lower panels) 
also displays a number of interesting features. First, the 
RWI was excited in the inner disk as early as 5 orbits. In 
contrast, the control run without solids (upper panels) has 
the instability appearing first in the outer disk, and at later 
times (10 orbits). The conclusion is that the particles in- 
duce vorticity on their own. Even though it is clear that 
this behavior has to do with free energy being transfered 
from the particle motion to the gas motion, it is not ob- 
vious if this result can be linked to the streaming insta- 
bility (Youdin & Goodman 2005) since the solids-to-gas 
ratio is not nearly as high as the one needed to excite it 
(pp/pg > 1). Instead, it is more likely that the backreaction 
is modifying the dispersion relation of the RWI. 

Another interesting feature of this run is that although 
the RWI was excited in the inner edge of the dead zone 
as early as 5 orbits, the outer edge just went unstable as 
late as 15 orbits. In contrast, the control run (upper panels) 
shows the outer disk going unstable at 10 orbits. As the 
backreaction hastens the growth of the RWI in the inner 
edge, it is unclear why it should stall it in the outer edge. 
One possibility is that the Rossby waves launched by the 
edge that first goes unstable interferes destructively with 
the perturbations fighting to grow in the other edge. 



14 



Lyra et al.: Planet formation bursts 



100 



Maximum drift (xQ=1) - Z =300 g cm" 3 - q=0.5 



o 40 




100 



80 



60 



Maximum drift (x£2=1 ) - L =1 50 g cm 3 - q=1 .5 



General 
Epstein 
Stokes 

mean free path 




Fig. 12. The radius of the particle subject to maximum drift (tQ=1) for our choice of parameters (left panel) and in the MMSN (right 
panel). q=—d In E/d lnr is the power law of the surface density profile. The profile is very flat compared to the ones predicted by the 
limiting cases of Epstein and Stokes drag, especially for our choice of parameters. The vortices in the inner and outer edge of the 
dead zone should have similar efficiency on trapping particles of a given size. 



The dominant mode also changed from the dusty to 
the dustless case. The latter has ra=4 and m=5 modes be- 
ing dominant in the inner edge. In the dusty case a m=2 
mode is seen instead. However, since in the dustless case 
it is the outer edge that displays am=2 mode, another ex- 
planation comes to mind. As the models are 2D, we are 
probably witnessing the inverse cascade phenomenon due 
to enstrophy conservation. The vortices are simply cascad- 
ing energy towards the larger scales, so the edge that goes 
unstable first (outer in the dustless case, inner in the dusty) 
will also reach a dominant m=2 mode first (possibly also 
m=l at later times). 



If this is the case, then self-gravity somehow halts 
the inverse cascade that took place in the dustless non- 
selfgravitating case. It is also instructive to compare the 
dusty non-self gravitating run (lower panel) with the dusty 
self gravitating run of LJKP08 (Fig. 1 of that paper). In that 
case, the ra=4 was dominant in the outer edge of the dead 
zone until the end of the simulation at 200 orbits. We also 
perform a test (Fig. 15) that consists of switching the self- 
gravity off in the run of LJKP08, and checking the evo- 
lution of the vortices. Without self -gravity, the ra=4 mode 
turns into a ra=3 mode in less than 15 orbits. In the inner 
edge, a m=2 mode developed out of the otherwise domi- 
nant m=3 mode. The inverse cascade indeed resumed. 



This result was also very recently reported by 
Mamatsashvili & Rice (2009). Without self-gravity, the vor- 
tex size is limited by the pressure scale height H. Once 
vortices grow to sizes of a few times H, the vortical flow 
becomes super-sonic. The vortex then radiates density 
waves that carry energy away and limit further growth. 
Mamatsashvili & Rice (2009) point that in the presence 
of self-gravity, the Jeans length Aj ~ QH, where Q = 
(kc s )/(GtcE) is the Toomre Q parameter, poses another 
limitation to the maximum size of a vortex. In Fig. 15 
we see that the vortices in the self-gravitating runs show 
Q « 1. The growth seen when self-gravity is switched off 
is a result of this constrain being lifted. 



7. Limitations of the model 

The presented models are admittedly simplified. In this 
section, we state what we consider the main limitations 
of our calculations to be. 

7. 1 . Two-dimensionality 

The most stringent limitation of the models is the 2D ap- 
proximation, which leads to a number of features, stated 
below. 

7.1.1. Vortex formation and survival 

The question of the excitation and sustainability of vor- 
tices in three dimensions is the matter of an old, yet un- 
settled, debate. Once excited, anticyclonic vortices are eas- 
ily maintained in 2D simulations where, unless viscosity is 
present, they cannot decay and will instead merge, grow- 
ing in size in a cascade of energy towards the largest scale 
of box (e.g., Johnson & Gammie 2005). However, three- 
dimensional studies in the context of protoplanetary disks 
found that tall vortex columns are destroyed, both in non- 
stratified (Shen et al. 2006) and in stratified (Barranco & 
Marcus 2005) local boxes. This phenomenon is understood 
as a result of the elliptic instability (Crow 1970, Gledzel 
et al. 1975, Kerswell 2002), by which the stretching term 
(cv • V)n, absent in 2D, breaks down elliptical stream- 
lines such as vortical flow. For a vortex to grow in 3D, the 
baroclinic term p~ 2 W p x Vp has to counter the stretching 
term. 

An indication that vortices can be sustained in three di- 
mensions is present in the study of Edgar & Quillen (2008). 
These authors simulate a stratified disk in spherical co- 
ordinates with an embedded giant planet. In their invis- 
cid run, the RWI is excited, leading to Rossby vortices at 
the edges of the gap, much like as in the 2D runs of de 
Val-Borro et al (2007). The vortices launched in three di- 
mensions are long lived and vertically extended, appar- 
ently following the same scale height as the surround- 
ing disk. We remark that the MRI-generated vortex of 
Fromang & Nelson (2005) is also seen to be long-lived in a 
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Fig. 13. Time series of the gas and solid density (for each indi- 
vidual particle size), showing the destruction of the embryos of 
10 cm particles when the inner vortices undergo the transition 
from the m=3 to the m=2 mode. At 75 orbits (upper panels), em- 
bryos are seen in both the inner and outer vortices, for both 10 cm 
and 30 cm particles. At 85 orbits one of the embryos in the 10 cm 
phase was disrupted, followed by a second at 95 orbits, and the 
last one ten orbits later. The embryos composed of 30 cm particles 
also experience tides. But as their gravitational Weber numbers 
are smaller, they just undergo splitting, the large fragments be- 
ing more stable against erosion than the embryos composed of 
particles of 10 cm. 



unstratified global disk. The studies of Edgar & Quillen 
(2008) and Fromang & Nelson (2005) both use a locally 
isothermal equation of state, which has large-scale non- 
zero baroclinity due to the static radial temperature gradi- 
ent. Furthermore, the existence of the RWI in 3D is demon- 
strated by the simulations of Meheut et al. (2008). 



7.1.2. Strength of the vortices 

A major impact of the 2D approximation is the inverse 
cascade due to enstrophy conservation that overpowers 
the vortices. An in depth study of the formation, devel- 
opment and structure of Rossby vortices in 3D global ac- 
cretion disks is needed to realistically address the issue of 
planet formation inside these structures. 



Fig. 14. Upper panels. Evolution of a disk without particles and 
without self-gravity, which serves as a control run for the next 
plots. 

Middle panels. Evolution of a disk without particles but with self- 
gravity. The difference compared to the upper panels is that the 
dominant mode in the outer disk changed from m=2 to m=5. Self- 
gravity modifies the dispersion relation of the RWI, or it stalls 
the inverse cascade of power known to occur in 2D turbulence, 
or both. 

Lower panels. Evolution of a disk without self-gravity but with 
particles. The backreaction leads to an early excitation of the RWI 
in the inner edge of the dead zone. Conversely, the outer edge 
goes unstable later when compared to the other two runs. Since 
the particle density is not high enough to excite the streaming 
instability, we take it as evidence that the backreaction modifies 
the dispersion relation of the RWI. 



7.1.3. Particle sedimentation 

Another limitation posed by the two-dimensionally is that 
the particles and the gas have the same infinitely thin scale 
height. The result of this is that the back-reaction of the 
drag force from the particles onto the gas is largely un- 
derestimated in our models. In 3D disks, the midplane 
particle layer is far denser due to sedimentation, so the 
ratio Pp/pg is far greater than the ratio E p /E g used in 
Eq. (3). The stronger backreaction that ensues is known 
to excite the streaming instabilities if pp/p g > 1 (Youdin 
& Goodman 2005, Youdin & Johansen 2007, Johansen & 
Youdin 2007). This instability enhances particle clumping, 
thus aiding collapse (Johansen et al. 2007). However, the 
effect of this strong backreaction on the vortices is poorly 
known. 



7.1.4. Different particle scale heights 

As the particles sediment, what sets the particle scale 
height is the equilibrium between turbulent diffusion and 
vertical gravity. Controlled by the drag force, the turbu- 
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lent diffusion depends on the particle radius, and so does 
the equilibrium scale height of the solids (Dubrulle et al. 
1995). Because of this, particles of lm radius settle in a 
thinner layer than those of 1cm particles. Inside a vor- 
tex, turbulent motions are expected to be weaker (Klahr 
& Bodenheimer 2006), bringing the layer of solids closer 
to a 2D configuration, but a dependence on radius is still 
expected. We could not model this effect on our simula- 
tions. 

7.1.5. Gas tides and mass loss 

The strength of the disrupting gas tides is yet another ef- 
fect related to the difference between 2D and 3D mod- 
els. As discussed in Sect. 3.2, the tides are proportional 
to the gas-to-solids ratio pg/ ' p v , thus expected to be much 
weaker in 3D where sedimentation considerably increases 
Pp in the midplane. As the vorticity is also expected to be 
weaker in 3D, the peak of p g at the vortex's eye would be 
weaker than in a 2D calculation, further weakening the ef- 
fect of tides. We therefore anticipate embryos formed in 3D 
simulations to be significantly less prone to mass loss than 
the ones presented in this study. 

7.2. Equation of state 

In this study, we used very simple equations of state: 
isothermal (Eq. (7)) or adiabatic (Eq. (24)). The effect of 
the equation of state can be appreciated by seeing the 
evolution of the Solberg-Hoiland criterion in isothermal 
and adiabatic simulations. In the former, it is the epicyclic 
frequency that brings k 2 + N 2 to negative values, while 
in the latter the criterion is broken mainly by the Brunt- 
Vaisala frequency. The excitation of the Rossby wave insta- 
bility is greatly favored in the presence of a strong entropy 
gradient, and made more difficult (yet not impossible) as 
the disk approaches isothermality. Therefore, an accurate 
modeling of the energy budget - solving for radiative cool- 
ing and turbulent heating -, is something to pursue in or- 
der to more realistically address the evolution of the RWI 
and the issue of planet formation inside the vortices that 
constitute its saturated state. 

7.3. Aerodynamics of the embryo 

The aerodynamics of a super-particle is controlled by the 
radius a. of the individual rocks. This means that even 
though the ensemble of rocks has the same position and 
velocity, there is still space between them and therefore 
they have contact with the gaseous nebula through all 
their surface area. 

However, after gravitational collapse occurs, the solids 
are not any longer an ensemble of pebbles and boulders 
with free space between them, but a single massive object 
of planetary dimensions. This leads to a radical change in 
the aerodynamical properties. Yet, in our simulations, we 
still consider the collapsed body as an ensemble of super- 
particles, with the aerodynamical properties of individual 
pebbles and boulders. This is certainly a limitation of the 
model. 

For a large planet, the correct treatment would be 
to consider that, after collapse, we leave the regime of 
particle-gas Stokes drag and enter planet-disk interaction 



by gravitational friction (type I migration). In the solar 
nebula the two drags have similar strength for bodies of 
100 km. As we solve for the self-gravity of the gas, the lat- 
ter is included in our model, albeit limited by the resolu- 
tion of the grid. The fact that we keep using Epstein-Stokes 
drag on the super-particles after collapse might make an 
embryo more stable, especially in view of the very effec- 
tive dynamical cooling provided by the drag force (Fig. 3). 
Substituting collapsed clusters by a sink particle that feels 
the gas gravity but not the gas drag is a possible solution, 
but also has caveats on its own. The evolution of sink par- 
ticles depends too much on artificial numerical parameters 
such as the accretion radius. Furthermore, a sink particle 
does not suffer tidal effects, which we showed to be non- 
negligible. 

7.4. Coagulation and Fragmentation 

As dust particles are drawn together, electromagnetic in- 
teractions occur at their surface, causing sticking under fa- 
vorable conditions. Brauer et al. (2008b) show that den- 
sity enhancements like the ones we see - where matter 
accumulates due to a discontinuity in viscosity -, dra- 
matically favor coagulation. As particles are drawn to- 
gether and the relative velocities are reduced, growth by 
coagulation occurs for a range of mass accretion rate m 
and threshold fragmentation velocity %. They find that 
the meter size barrier can be breached for mass accretion 
rates up to ra=10 _8 Mo/yr (for %=10ms _1 ) and thresh- 
old fragmentation velocities no less than i;f t =5ms -1 (for 
m=S x lO _9 M /yr). 

This raises the possibility that even before the RWI ex- 
cites the vortices, coagulation will have depleted the pop- 
ulation of centimeter and meter sized objects onto bodies 
that are too large for our proposed mechanism to occur 
efficiently. As we see, it is preferentially the 10 and 30 cm 
sized particles that concentrate into planetary embryos. 

The timescale for coagulation, however, is much longer 
than the time-scale for gravitational collapse. We see 
growth to Mars size taking place in only five orbits (~60 
yr), while growth by coagulation from meter to kilometer 
size occurs on timescales of a few thousand years accord- 
ing to Brauer et al. (2008b). On the other hand, it could 
as well be that the favorable environment provided by the 
vortices act as to speed up coagulation even further. This, 
of course, is not bad. Growth beyond the preferred size 
will lead to decoupling from the gas and ejection from the 
vortex that, in the end, behaves as a planetesimal factory. 

Fragmentation is an important piece of reality that we 
did not include in our model. Nevertheless, we showed 
in Sect 4.5 that the majority of the particles were never 
involved in collisions with speeds greater than 1ms -1 . 
These are of course very good news for planet formation. 
However, we feel the need to stress that the time inter- 
val between snapshots in Fig. 9a is of whole orbits (to- 
taling 200 snapshots). The number of high-speed impacts 
could be greater had we checked the collision speeds at 
every time-step instead. Although desirable, this would 
have made the simulations computationally very expen- 
sive since it must be done in runtime. The result of Fig. 9b 
should therefore be taken only as further evidence for low 
collisional speeds inside vortices, not as conclusive proof 
of it. Carballido et al. (2008) further point that the low colli- 
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sional speeds at low separation may be unrealistic. This is 
because the particles couple to the smallest eddies, whose 
size is a function of the mesh Reynolds number. These au- 
thors find average collisional speeds of 0.05c s for particles 
of stopping time 0^1^=0.2 (a. ^15 cm in our models), but 
notice a sharp decrease of the collisional speeds towards 
smaller separations. In our simulations, we are consider- 
ing encounters that happen inside a grid cell, where we do 
not resolve the velocity field, so this may indeed be the rea- 
son behind the low collisional velocities we find. However, 
we point that there is a difference between the simulations 
of Carballido et al. (2008) and those presented in this pa- 
per. They considered particle concentrations in the tran- 
sient pressure maxima of the turbulence generated by the 
MRI, whereas we consider particle concentrations within 
long lived anticyclonic vortices. As vortex structures tend 
to merge and grow, there is less power available at the 
smallest scales when compared to MRI turbulence. 



7.5. Disk mass 



We stated in Sect 2.2 that the constrain of having 30 Earth 
masses of solid in the inner disk, together with the inter- 
stellar dust-to-gas ratio of e=10 -2 , and a surface density 
slope of -0.5 implies 0.01 solar mass within the modeled 
range. This disk can be considered massive since the shal- 
low slope then implies a massive outer disk: 0.3 solar mass 
within 100 AU. The actual mass is somewhat lower be- 
cause the disk goes isothermal after a given radius. If the 
temperature at 5.2 AU is 100K, it will have fallen to 10 K at 
50 AU and below the cosmic microwave background tem- 
perature at 200 AU. As the surrounding temperatures of 
the molecular cloud are of ^10-20 K, the 1/r fall of the 
temperature has to stop at some point. Boss (2002) and 
Pinotti et al. (2005) account for this by using a tempera- 
ture profile of T = T r _ ^ T + 7g, where the parameter 7g 
works as a background temperature (see also Papaloizou 
& Terquem 1999), leading to an essentially isothermal pro- 
file after a certain distance. The constancy of the mass ac- 
cretion rate then dictates that the slope of the surface den- 
sity has to change accordingly, to -1.5. This, in turn, leads 
to a lower disk mass. Taking into account this isothermal 
profile in the very outer disk, using Tg=10K, and the cor- 
responding change of slope in Z, the disk sets at a constant 
value of Q=l after 50 AU, with a mass of 0.25 solar mass 
at 100 AU. If T B =20K, the mass decreases to 0.16 M© at 
100 AU, with Q=3 after 25 AU. 

Observational studies in contrast obtain a distribution 
of disk masses with median of 5 x 10 -3 M© (Andrews & 
Williams 2005). However, considering an interstellar dust- 
to-gas ratio of 10 -2 , such disks contain only 15 M© of 
solid material, and cannot form the solar system. We take 
it as an indication that or the masses are underestimated 
(Hartmann 2008) or that, as these observations correspond 
to relatively older disks (1-4 Myr), the gas has already dis- 
sipated significantly (Hueso & Guillot 2005). Our model 
then should be representative of young disks, therefore 
corresponding to an early formation of relatively massive 
planetary embryos. 



8. Summary and conclusions 

We have undertaken simulations of disks of gas and solids, 
where the solids are represented by Lagrangian particles 
of radii 1,10, 30, and 100 cm. We show that in the borders 
of the dead zone, where the accretion flow stalls due a dif- 
ference in turbulent viscosity, rapid gravitational collapse 
of the particles into protoplanets occurs within the vortices 
that form due to the excitation of the Rossby wave insta- 
bility (Lovelace et al. 1999). As shown in LJKP08, over 300 
gravitationally bound planetary embryos were formed, 20 
of them being more massive than Mars. The mass spec- 
trum follows a power law of index — 2.3±0.1 in the inter- 
val -2.0<log(M/M©)< -1.2. 

Although for the main results of this study we have 
used sharp viscosity jumps to model the transition be- 
tween the active and dead zones, we show that the RWI is 
excited up to viscosity jumps as smooth as Ar=2H where 
H is the pressure scale height. For this conclusion, we used 
the Solberg-Hoiland criterion as providing a conservative 
estimate of whether the RWI is excited. The consequence 
of increasing the width of the viscosity jumps seems to be 
that the threshold of instability is reached at increasingly 
longer times. It only takes five orbits with Ar=0.2H (the 
usual width used in the models presented in this paper), 
but takes 350 orbits for Ar=2H. We also assessed if the vor- 
tices would survive in the more realistic environment of a 
turbulent disk by making the location of the viscosity shift 
oscillate with an amplitude typical of the scale length of 
MRI turbulence over the period of a Keplerian revolution. 
We find that this has little effect on the excitation of the 
RWI and saturation into vortices. 

We model the solid phase with Lagrangian superpar- 
ticles representing physical pebbles and rocks of different 
size (1, 10, 30, and 100 cm). As these particles are subject to 
different gas drag, an aerodynamical sorting by size takes 
place within the vortices. The consequence of this is that 
the first bound structures are formed of single particles 
species. This is a very interesting result, since it is an ob- 
served fact that planetesimals are formed of similar-sized 
building blocks (Scott & Krot 2005). These building blocks 
seem to be sub-mm sized grains, but different nebula pa- 
rameters could work as to trap smaller particles. Youdin 
(2008) also points that the stationary point of a particle 
trapped in vortical motion is shifted azimuthally with re- 
spect to the eye, according to its radius a.. We indeed 
see that clumps of particles of different size, which col- 
lapse into different embryos inside the same vortex, usu- 
ally retain significant azimuthal displacements between 
each other for long times instead of forming a single, more 
massive, embryo at the vortex eye. This may or may not 
be a result of the size-dependent azimuthally shifted sta- 
tionary points of Youdin (2008). 

A collapsed embryo is observed to be very compact. 
The compactness is mainly provided by the drag force, 
which makes the system very dissipative (dynamical cool- 
ing by gravity alone works on much longer timescales). 
Collapse towards // infinite ,/ density is seen to occur in 
some cases, with most of the particles occupying the same 
position in space (limited by numerical precision). In the 
specific case when the particles dominate the time-step, 
the Courant condition leads a particle to overshoot the 
center of the mass, so that it executes oscillations about it, 
which in turn leads to a finite rms radius. We also observe 
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Fig. 15. In the simulation shown in LJKP08, the vortices maintained a ra=4 mode in the outer edge of the dead zone until the end of 
the simulation. By switching self-gravity off, we see that in less than 15 orbits the outer vortices merged into a m-3 mode, while the 
m=3 mode in the inner disk edge into a m=2 mode. This is because as Q decreases, the size of the vortices approaches the Jeans length 
scale, which effectively halts the inverse cascade of energy. 



that a clump of particles is susceptible to the disrupting ef- 
fects of ram pressure erosion and gravitational tides from 
the gas. Both effects are proportional to the local gas-to- 
solids density ratio. When the vortices in the inner border 
of the dead zone undergo the transition from the m=3 to 
the m=2 mode, their surface density increases, with dras- 
tic consequences for the embryos within them. The most 
massive embryo by that time, a protoplanet 6.7 times the 
mass of Mars, mostly formed of a. =30 cm particles, was 
split into two smaller objects, of 5.9 and 0.8 MMars/ due to 
the action of the gas tides. The fate of the embryos formed 
mostly of 10 cm was more dramatic. As the 10 cm particles 
experience stronger drag forces, the ram pressure is also 
stronger. During the mode transition, the combined effects 
of tides and erosion completely obliterated these embryos, 
leaving extended arcs of particles that did not reaccumu- 
late until the end of the simulation. We anticipate that this 
effect will be very reduced in 3D simulations. In 2D simu- 
lations, the ratio of the vertically integrated solids density 
to the gas column density Ep / Eg never gets much above 
unity even for the most massive embryo. In contrast, the 
ratio of the bulk density of solids to the volume gas den- 
sity pp/pg is greatly increased in the midplane of 3D disks 
due to sedimentation. 

We also observe that the solids modify the evolution 
of the RWI. We are drawn to this conclusion because a 
simulation without the backreaction of the drag force from 
the particles onto the gas developed vortices at later times 
when compared to the ones that include particle feedback. 
We stress that this is not due to the streaming instabil- 
ity, since the solids-to-gas ratio was much lower than the 
value needed to excite it (pp/pg > 1). Instead, it is more 
likely that the backreaction of the drag force contributes 



non-negligibly to the dispersion relation of the RWI. Self- 
gravity is also seen to play a role on modifying the evo- 
lution of the turbulence. We observe that in simulations 
without self-gravity, the disk tends to show less vortices 
at later times. In a simulation where we switched off the 
self-gravity after the vortices had developed, the domi- 
nant ra=4 mode in the outer edge of the dead zone rapidly 
turned into a m=3 mode. The vortices in the inner edge 
also quickly turned from a dominant m=3 mode to dis- 
playing a m=2 mode instead. We measured the Toomre Q 
parameter and found that the vortices have Q ~ 1. This 
constitutes further evidence that in the presence of self- 
gravity, vortex growth is not only limited by the pressure 
scale height but also by the Jeans length (Mamatsashvili & 
Rice 2009). 

An important finding in this paper is that under vor- 
tex trapping, the collisional speeds between particles are 
greatly reduced. We measured the collisional velocity his- 
tory of every particle that is bound to the most massive 
embryo at the end of the simulation, and found that the 
vast majority of them never experienced close encounters 
at speeds greater than 1ms -1 . This is well below the frag- 
mentation threshold and lends support to the long-held 
idea that vortices provide a superbly favorable environ- 
ment for planetary growth. Growth by coagulation be- 
yond the optimal size for planet formation is also avoided 
because the timescales for coagulation are much longer 
than the rapid timescale for gravitational collapse wit- 
nessed in our models. This does not exclude the possi- 
bility that coagulation itself is sped up within a vortex. 
In this case, the vortex will behave as a planetesimal fac- 
tory, quickly producing kilometer sized bodies that leave 
the vortex due to their decoupling from the gas. This, as 
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noted by Klahr & Bodenheimer (2006) is very positive 
for planet formation, since the formed planetesimals are 
then scattered through the disk, where they can be used 
to form planets independently of a vortex. Even though 
it implies that we are facing the comfortable position of 
a win-win situation for planet formation, one has to de- 
cide which process (planet formation by direct gravita- 
tional collapse or planetesimal formation by fast coagu- 
lation) is getting the upper hand inside the Rossby vor- 
tices. A definite answer to this question can only be drawn 
from a simulation that includes the processes of coagula- 
tion/fragmentation. Unfortunately, inclusion of sophisti- 
cated coagulation/ fragmentation models such as that of 
Brauer et al. (2007) would render a hydrodynamical sim- 
ulation overly expensive. A possible alternative would be 
a Monte Carlo model of dust coagulation, such as the one 
recently developed by Ormel et al. (2007). A further devel- 
opment of the Monte Carlo technique is described in Zsom 
& Dullemond (2008). The main difference between the two 
models is that while Ormel et al. (2007) simulate coag- 
ulation between real dust particles, Zsom & Dullemond 
(2008) use superparticles to model coagulation and frag- 
mentation. Therefore the latter one is more suitable for hy- 
drodynamical simulations such as ours and simple esti- 
mations show that this model could be adapted to a hydro 
model with no prohibitive overhead. 

We reiterate that the models presented suffer from a 
number of limitations, detailed in Sect 7. Some of them, 
like refining the particle mass resolution to the individ- 
ual pebbles and rocks, are beyond the capabilities of the 
current generation of computer models. Others, however, 
such as inclusion of detailed thermal physics, could be 
tackled with relatively little effort. We urge researchers ac- 
tive on the field to consider these problems. It is our hope 
that a coherent picture of planet formation in the magneti- 
cally dead zones of accretion disks shall emerge as a result 
of it. 
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Appendix A: Resolution study 

In Fig. A.l we present a resolution study of our models. 
The upper panels show Cartesian models, including self- 
gravity and solids, while the lower panels show cylindri- 
cal models with dustless non-self gravitating gas. Both are 
shown in selected snapshots at 10, 20, and 30 orbits. 

In the Cartesian runs, the vortices in the high- 
resolution run (512 2 ) are excited earlier than in the mid- 
dle resolution run (256 2 ). At later times, the vortices in the 
high-resolution run also appear sharper. While in the two 
runs the vortices in the outer edge of the dead zone look 
remarkably similar, the inner edge displays very different 
behavior. Up to 30 orbits the middle resolution run has not 
shown signs of prominent vortices. In contrast, the high 
resolution run had the inner edge developing vortices as 
early as 10 orbits. This is not only due to the high resolu- 
tion run having a wider inertial range, but also because a 
flow with cylindrical symmetry is more coarsely resolved 
near the origin in a Cartesian grid. Because of this, the 
most unstable modes are under-resolved in the inner disk 
of the middle resolution run. A Cartesian run with low res- 
olution (128 2 , not shown) did not develop vortices even in 
the outer disk by the same time. At 30 orbits, the density 
inside the vortices peak at similar values, E/Eq=33 and 
3.7 for the middle and high resolution runs, respectively. 

The cylindrical runs also illustrate the small amount 
of differences between the vortices in different runs. With 
better azimuthal resolution, the vortices are excited even 
in the 128 2 run, and in both runs they peak with surface 
density E/Eq=4.6. In fact, the main effect of resolution ap- 
pears to be that, as it increases, the RWI is excited increas- 
ingly earlier. The high-resolution run displays weaker vor- 
tices than the others because in this case we were forced to 
use shock viscosity. 



Appendix B: Emulating turbulent motions 

In this study, we considered the dead zone to be repre- 
sented by a static viscosity profile. In a more realistic sce- 
nario, turbulent motions caused by the MRI and variations 
in the coupling between the magnetic field and the plasma 
will give rise to a turbulent resistivity. This is expected to 
cause the border of the dead zone to vary in space and 
time, with implications for the evolution of the RWI. 

To assess the impact of space and time variability of 
the edges of the dead zone, we model the viscous jumps 
using Ar=0.01 but make the jumps shift radially in time by 
substituting r\ and r<i in Eq. (13) by 



log 10 (z/I (l _ 0) ) 



n(t) n [i + h sin (n K (r x )t)] 

r 2 (t) r 2 [1 + h sin (Ct K (r 2 )t)] 



(B.l) 
(B.2) 



where h=H/r is the aspect ratio. So, the location shifts by 
two scale heights over a Keplerian revolution. The results 
are shown in Fig. B.l, where we show the appearance of 
the disk at 30 orbits and the azimuthal average of (k 2 + 
N 2 )/Q 2 K (in the same 2D model, as opposed to ID as in 
Sect. 3.2). 

In this simulation, the RWI is still excited and vor- 
tices are launched. The main difference when compared 
to simulations with static profiles is that the instability 
takes more time to violate the Solberg-Hoiland criterion, 
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Fig. A.l. The development of the dead zone vortices in different 
resolution and grid geometry. With increased resolution the RWI 
is excited increasingly earlier. In general the vortices in cylindri- 
cal runs look sharper than in the Cartesian ones, due to the bet- 
ter azimuthal resolution. The Cartesian run of middle resolution 
(256 2 ) has too coarse azimuthal resolution in the inner disk, and 
only developed vortices in the inner edge of the dead zone at 
later times 40 orbits). Apart from these differences, the vor- 
tices look remarkably similar, having nearly the same peak den- 
sity and same vorticity. 



« 10 orbits, compared to 5 in the static case. This is due 
to the fact that the shifting viscosity jump smears the pres- 
sure maximum, so the amplitude of the pressure jump is 
shorter and the width is larger than in the static case. 



Appendix C: Onset of erosion 

According to Eq. (30), a tight distribution of particles un- 
der Epstein drag should suffer erosion if 
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> 1 
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Fig.B.l. Evolution of the RWI with a time-varying viscosity jump as specified by Eq. (13) and Eqs. (B.1)-(B.2). The panel on the 
left-hand side measures the violation of the Solberg-Hoiland criterion. The right-hand side panel shows the appearance of the disk 
at 30 orbits. Both the inner and outer edge quickly reach the threshold of instability (left panel). At 30 orbits, the inner edge already 
reached a saturated state and launched vortices (right panel). 
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Fig. C.l. A clump of 1000 massive particles moving against headwinds of Ma=0.5 and Ma=0.05, for different values of the friction 
time tQ. In the case of Ma=0.5, self-gravity cannot hold the clump together for rO < 10~ 2 . In our simulations it corresponds to 1 cm 
sized particles, approximately. Clumps formed of larger particles to not experience erosion. For the more subsonic motion, the effect 
of ram pressure is lower so the clump of rn=10 -2 is more stable against erosion. The case of rn=1.0 takes longer to contract because 
of the weaker drag force, which provides less dynamical cooling than in the case with rn=0.1. It eventually shrinks, as seen in the 
time series (bottom plot). 
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In this appendix, we perform numerical simulations 
to test the validity of this condition. We model a clump 
of 10 3 particles suffering a strong headwind from the gas 
(Ma=l/2). The blob of particles is initially set as a tight 
Gaussian distribution about the center of the grid, with 
surface density peak of Zp/Z$ = 2.31. In units where 
G=c s =Zq=Qq=1, the initial gravitational Weber number at 
the surface of the blob is therefore Wee - 7x 10~ 2 / (rO). 

We plot in Fig. C.l the evolution of clumps for four 
different values of rO. For t\Q=1.0 and r 0=0.1, Wee is 
less than 1 so the clump is stable against ram pressure 
and contracts. The other clumps (rf^O.Ol and r.O^O.001) 
have Wee above unit, and experience intense erosion. We 
also considered a flow of Mach number Ma=0.05. In this 
case the initial gravitational Weber number is Wee — 
7 x 10 -3 / (tO), and the clumps of smaller particles are 
supposed to be more stable. Indeed, this is what we see 
in the figure. The clump of particles of rf^O.Ol is now 
marginally stable and contracts, experiencing much less 
erosion than in the Ma=0.5 case. 

We would like to draw attention to an interesting fea- 
ture of Fig. C.l. The cases of rO=l and rf^O.l (upper 
panels) provide yet another perspective for the action of 
drag force cooling. Both clumps are stable against erosion 
but the clump of rf^O.l has shrunk considerably more 
than the clump of t 0=1.0, which looks very extended. 
What is happening is that the clump with rO=l is too 
weakly coupled to the gas and therefore takes longer to 
collapse, as seen in the time series of the rms position 
(Fig. C.l, lower panels). 



